/*
 heat-spectral lib. solves time-dependent PDEs using SEM-NI method.
 Copyright (C) 2009, Altieri Michele and Stefanoni Andrea
 Copyright (C) 2008, Proverbio Alessandro and Radice David

 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.
 
 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 GNU General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

/*!
    \file umfmatrix.cpp
    \brief This file contains the UMFMatrix class

*/

#include <cassert>
#include <iostream>

#include "umfmatrix.hpp"

namespace spectral
{

UMFMatrix::UMFMatrix( void ) : _Ap(NULL),
                               _Ai(NULL),
                               _Ax(NULL),
                               _sizeA(0),
                               _sizeNZ(0),
                               _Symbolic(NULL),
                               _Numeric(NULL),
                               _Control(NULL),
                               _Info(NULL),
                               _sys(UMFPACK_Aat),
                               _isFactorized(false) {
    _Control=new double[UMFPACK_CONTROL];
    umfpack_di_defaults(_Control);
    _Info=new double[UMFPACK_INFO];
}



UMFMatrix::UMFMatrix( int sizeA, int sizeNZ ) : _sizeA(sizeA),
                                                _sizeNZ(sizeNZ),
                                                _Symbolic(NULL),
                                                _Numeric(NULL),
                                                _Control(NULL),
                                                _Info(NULL),
                                                _sys(UMFPACK_Aat),
                                                _isFactorized(false) {
    assert(_sizeA>0);
    assert(_sizeNZ>0);
    _Ap = new int[_sizeA+1];
    _Ai = new int[_sizeNZ];
    _Ax = new double[_sizeNZ];
    _Control=new double[UMFPACK_CONTROL];
    umfpack_di_defaults(_Control);
    _Info=new double[UMFPACK_INFO];
}

UMFMatrix::UMFMatrix( UMFMatrix const & oldUMFMatrix ) : _sizeA(oldUMFMatrix._sizeA),
                                                         _sizeNZ(oldUMFMatrix._sizeNZ),
                                                         _Symbolic(NULL),
                                                         _Numeric(NULL),
                                                         _Control(NULL),
                                                         _Info(NULL),
                                                         _sys(UMFPACK_Aat),
                                                         _isFactorized(false) {

    assert(_sizeA>0);
    assert(_sizeNZ>0);

    _Ap = new int[_sizeA+1];
    _Ai = new int[_sizeNZ];
    _Ax = new double[_sizeNZ];

    for(int i(0); i<=_sizeA; ++i )  {
        _Ap[i] = oldUMFMatrix._Ap[i];
    }

    for(int i(0); i< _sizeNZ; ++i )  {
        _Ai[i] = oldUMFMatrix._Ai[i];
        _Ax[i] = oldUMFMatrix._Ax[i];
    }
    _Control=new double[UMFPACK_CONTROL];
    umfpack_di_defaults(_Control);
    _Info=new double[UMFPACK_INFO];
}

UMFMatrix & UMFMatrix::operator=(UMFMatrix & B) {
    delete [] _Ap;
    delete [] _Ai;
    delete [] _Ax;

    _sizeA = B.getsizeA();
    _sizeNZ = B.getsizeNZ();

    _Ap = new int[_sizeA+1];
    _Ai = new int[_sizeNZ];
    _Ax = new double[_sizeNZ];

    for(int i(0);i<=_sizeA;++i) {
	_Ap[i] = B.getAp(i);
    }

    for(int i(0);i<_sizeNZ;++i) {
	_Ai[i] = B.getAi(i);
	_Ax[i] = B.getAx(i);
    }
    return *this;
}

UMFMatrix::~UMFMatrix()  {
    clean();
    delete [] _Ap;
    delete [] _Ai;
    delete [] _Ax;
    delete [] _Control;
    delete [] _Info;
}

void UMFMatrix::resize_sizeA(int sizeAnew) {
    _sizeA=sizeAnew;
    delete [] _Ap;
    _Ap=new int[_sizeA+1];
    if(_isFactorized) {
        clean();
    }
}

void UMFMatrix::resize_sizeNZ(int sizeNZnew) {
    _sizeNZ=sizeNZnew;
    delete [] _Ai;
    delete [] _Ax;
    _Ai = new int[_sizeNZ];
    _Ax = new double[_sizeNZ];
    for(int i=0;i<_sizeNZ;++i) {
        _Ax[i]=(double)0;
    }
    if(_isFactorized) {
        clean();
    }
}

void UMFMatrix::resize(int sizeAnew, int sizeNZnew)  {
    resize_sizeA(sizeAnew);
    resize_sizeNZ(sizeNZnew);
}

void UMFMatrix::solve(std::valarray<double> const & b,
                      std::valarray<double> & u) const {
    solve(&b[0],&u[0]);
}

void UMFMatrix::solve(double const * b, double * x) const {
    if(!_isFactorized) {
        std::cerr << "UMF: We can't solve if we don't factorize the matrix...\n";
        exit(1);
    }

    _status = umfpack_di_solve (_sys, _Ap, _Ai, _Ax, x, b, _Numeric, _Control, _Info) ;

    // Warnings
    if( _Info [UMFPACK_STATUS] != UMFPACK_OK ) {
        if(_Info [UMFPACK_STATUS] == UMFPACK_ERROR_invalid_Numeric_object )
            std::cerr << "UMF: Solve: Invalid numeric object " << std::endl;
        if(_Info [UMFPACK_STATUS] == UMFPACK_ERROR_out_of_memory )
            std::cerr << "UMF: Solve: Out of memory " << std::endl;
        if(_Info [UMFPACK_STATUS] == UMFPACK_ERROR_argument_missing )
            std::cerr << "UMF: Solve: Argument missing " << std::endl;
        if(_Info [UMFPACK_STATUS] == UMFPACK_ERROR_invalid_system )
            std::cerr << "UMF: Solve: Invalid system " << std::endl;
        else  std::cerr << "UMF: Solve is NOT ok " << std::endl;
    }

    if( _Info[UMFPACK_STATUS] != UMFPACK_OK ) {
        exit(1);
    }
}

int UMFMatrix::getsizeA() const {
    return _sizeA;
}

int UMFMatrix::getsizeNZ() const {
    return _sizeNZ;
}

UMFMatrix & UMFMatrix::operator+=(UMFMatrix & B) {
    assert(B.getsizeNZ() == _sizeNZ);
    for (int i=0; i<_sizeNZ; ++i) {
	_Ax[i] += B.getAx(i);
    }
    return *this;
}

UMFMatrix & UMFMatrix::operator*=(double const & alpha) {
    for (int i=0;i<_sizeNZ;++i) {
	_Ax[i] = alpha * _Ax[i];
    }
    return *this;
}

void UMFMatrix::printPattern(char const * name) const {
    std::fstream pattern;
    pattern.open(name, std::ios_base::out);

    int items = 0;
    int index = 0;
    int length = 0;

    for(int j=0;j<_sizeA;++j) {
	items = _Ap[j+1] - _Ap[j];
	for(int i=0;i<_sizeA;++i) {
	    if(index>=items) {
		for (int k=i;k<_sizeA;++k) {
		    pattern << "0 ";
		}
		break;
	    }
	    if(_Ai[length+index] == i) {
		pattern << "1 ";
		++index;
	    } else {
		pattern << "0 ";
	    }
	}
	length += index;
	index = 0;
	pattern << std::endl;
    }
}

void UMFMatrix::multiply(std::valarray<double> const & in,
			 std::valarray<double> & out) const {
    multiply(&in[0],&out[0]);
}

void UMFMatrix::multiply(double const * V, double * AV) const {
    for(int i=0;i<_sizeA;++i) {
	AV[i]=0.0;
    }
    for(int i=0;i<_sizeA;++i) {
	for(int j=_Ap[i];j<_Ap[i+1];++j) {
	    assert(j<_sizeNZ);
	    assert(_Ai[j]<_sizeA);

	    AV[i]+=_Ax[j]*V[_Ai[j]];
	}
    }
}

double * * UMFMatrix::fullmatrix() const {
    double * * A;
    A=new double * [_sizeA];
    for(int i=0;i<_sizeA;++i) {
        A[i]=new double [_sizeA];
    }

    for(int i=0;i<_sizeA;++i) {
        for(int j=0;j<_sizeA;++j) {
            A[i][j]=0.e0;
        }
    }

    int begin;
    int end;
    int row;
    for(int col=0;col<_sizeA;++col) {
        begin=_Ap[col];
        end=_Ap[col+1];
        for(int idx=begin;idx<end;++idx) {
            row=_Ai[idx];
            A[row][col]=_Ax[idx];
        }
    }
    return A;
}

void UMFMatrix::factorize() {
    // =================== Call symbolic ===================
    _status = umfpack_di_symbolic ( _sizeA, _sizeA, _Ap, _Ai, _Ax, &_Symbolic, _Control, _Info );

    // Troubleshooting
    if( _Info [UMFPACK_STATUS] != UMFPACK_OK ) {
        if(_Info [UMFPACK_STATUS] == UMFPACK_ERROR_internal_error )
            std::cerr << "UMF: Symbolic: A bug occours " << std::endl;
        if(_Info [UMFPACK_STATUS] == UMFPACK_ERROR_out_of_memory )
            std::cerr << "UMF: Symbolic: Out of memory " << std::endl;
        if(_Info [UMFPACK_STATUS] == UMFPACK_ERROR_argument_missing )
            std::cerr << "UMF: Symbolic: Argument missing " << std::endl;
        if(_Info [UMFPACK_STATUS] == UMFPACK_ERROR_invalid_matrix )
            std::cerr << "UMF: Symbolic: Error in the matix " << std::endl;
        else  std::cerr << "UMF: Symbolic is NOT ok " << std::endl;
    }

    // =================== Call numeric ===================
    _status = umfpack_di_numeric (_Ap, _Ai, _Ax, _Symbolic, &_Numeric, _Control, _Info) ;

    if( _Info [UMFPACK_STATUS] != UMFPACK_OK ) {
        if(_Info [UMFPACK_STATUS] == UMFPACK_WARNING_singular_matrix )
            std::cerr << "UMF: Numeric: Singular matrix " << std::endl;
        if(_Info [UMFPACK_STATUS] == UMFPACK_ERROR_out_of_memory )
            std::cerr << "UMF: Numeric: Out of memory " << std::endl;
        if(_Info [UMFPACK_STATUS] == UMFPACK_ERROR_argument_missing )
            std::cerr << "UMF: Numeric: Argument missing " << std::endl;
        if(_Info [UMFPACK_STATUS] == UMFPACK_ERROR_invalid_Symbolic_object )
            std::cerr << "UMF: Numeric: Invalid symbolic object " << std::endl;
        if(_Info [UMFPACK_STATUS] == UMFPACK_ERROR_different_pattern )
            std::cerr << "UMF: Numeric: Different pattern " << std::endl;
    }

    _isFactorized=true;
}

void UMFMatrix::clean() {
    umfpack_di_free_symbolic (&_Symbolic);
    umfpack_di_free_numeric (&_Numeric);
    _isFactorized=false;
}

UMFMatrix & operator+(UMFMatrix A, UMFMatrix B) {
    return A+=B;
}

UMFMatrix & operator*(UMFMatrix A, double alpha) {
    return A*=alpha;
}

UMFMatrix & operator*(double alpha, UMFMatrix A) {
    return A*=alpha;
}

} // namespace
